% This Matlab program generates figure 2;

clear
tic

load ('ForBootstraps');

reps = 50;

stcnty = double( st_dataset.stcnty );
n = size(stcnty,1);
S = sparse( [1:292308]', stcnty, ones(n,1) );  % (n x m) '

tram = double(st_dataset.tr_am);
availcredit = double(st_dataset.availcredit);
limit = double(st_dataset.limit);
pctblack = double(st_dataset.pctblack);
util_dol = double(st_dataset.util_dol);
bal = double(st_dataset.bal);
sc_black = double(st_dataset.sc_black);
unusedlim = double(st_dataset.unusedlim);

demographics = [ st_dataset.age st_dataset.gt_hs_male st_dataset.eq_hs_male st_dataset.gt_hs_female st_dataset.eq_hs_female st_dataset.nonmarried_male st_dataset.nonmarried_female st_dataset.widowed_male st_dataset.widowed_female st_dataset.divorced_male st_dataset.divorced_female st_dataset.foreignborn ];
income = [ st_dataset.income_10 st_dataset.income_15 st_dataset.income_20 st_dataset.income_25 st_dataset.income_30 st_dataset.income_35 st_dataset.income_40 st_dataset.income_45 st_dataset.income_50 st_dataset.income_75 st_dataset.income_100 st_dataset.income_125 st_dataset.income_150 st_dataset.income_200 ];
emp = [ st_dataset.publicassistance st_dataset.employment ];
housing = [ st_dataset.vacant st_dataset.ownocc st_dataset.withmort st_dataset.medianrent_b st_dataset.mhv ];

clear st_dataset;

n = size(availcredit,1);
[sort_fips,fips_index] = sort(stcnty);
chg_fips = zeros(n,1);
chg_fips(1) = 1;
chg_fips(2:n) = (sort_fips(2:n)~=sort_fips(1:n-1));
fips_grp = cumsum(chg_fips,1);
fips_grp(fips_index) = fips_grp;
fips_dummies = sparse(dummyvar(fips_grp));
%fips_dummies(:,1) = [];

% Baseline with interaction
y = limit;
x = [pctblack sc_black tram util_dol];

[n,k] = size(x);
x = [x fips_dummies];
k1 = size(x,2);

xpxi = (x'*x)\eye(k1);
beta = xpxi*(x'*y);  % the first two elements should be the PCT_BLACK variables
beta(1:2)

scorerange = linspace(min(tram),max(tram),100);
meanx = mean(x,1);
allwhite = meanx;
allwhite(1:3) = 0;
allwhite = (allwhite*beta) + ( beta(3) * scorerange);

allblack = meanx;
allblack(1:3) = [1 0 0];
allblack = (allblack*beta) + ( (beta(2) + beta(3))* scorerange );

figure(1);
clf;
subplot(1,2,1);
set(gcf,'PaperUnits','inches','PaperPosition',[0 0 6 3]);

plot(scorerange,allwhite,'k-',scorerange,allblack,'k--');
xlim([300,900]);
legend(['All White';'All Black'],'Location','SouthEast');
legend boxoff;
xlabel('Credit Score');
ylabel('Credit Limit (in thousands)');
ylim([0,60]);
title('A) Replication');
set(gca,'position',[0.1 0.14 0.35 0.7]);

% Baseline with interaction and income
y = limit;
x = [pctblack sc_black tram util_dol income];

[n,k] = size(x);
x = [x fips_dummies];
k1 = size(x,2);

xpxi = (x'*x)\eye(k1);
beta = xpxi*(x'*y);  % the first two elements should be the PCT_BLACK variables
beta(1:2)

scorerange = linspace(min(tram),max(tram),100);
meanx = mean(x,1);
allwhite = meanx;
allwhite(1:3) = 0;
allwhite = (allwhite*beta) + ( beta(3) * scorerange);

allblack = meanx;
allblack(1:3) = [1 0 0];
allblack = (allblack*beta) + ( (beta(2) + beta(3))* scorerange );

subplot(1,2,2);

plot(scorerange,allwhite,'k-',scorerange,allblack,'k--');
xlim([300,900]);
legend(['All White';'All Black'],'Location','SouthEast');
legend boxoff;
xlabel('Credit Score');
ylabel('Credit Limit (in thousands)');
ylim([0,60]);
title(['B) Replication without Distortion';'Effects Plus Income Controls     ']);
set(gca,'position',[0.6 0.14 .35 .7]);

% This will output the graphic image for figure 2;
print -dpng ScoreDifferenceGraphs.png -r100



